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Abstract 

The nucleation of cavities in a homogenous polymer under tensile strain is 
investigated in a coarse-grained molecular dynamics simulation. In order to 
establish a causal relation between local microstructure and the onset of cav- 
itation, a detailed analysis of some local properties is presented. In contrast 
to common assumptions, the nucleation of a cavity is neither correlated to a 
local loss of density nor, to the stress at the atomic scale and nor to the chain 
ends density in the undeformed state. Instead, a cavity in glassy polymers 
nucleates in regions that display a low bulk elastic modulus. This criterion 
allows one to predict the cavity position before the cavitation occurs. Even 
if the localization of a cavity is not directly predictable from the initial con- 
figuration, the elastically weak zones identified in the initial state emerge as 
favorite spots for cavity formation. 
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1 1. Introduction 

2 Under hydrostatic stress conditions, failure of amorphous polymers oc- 

3 curs through cavitation, often followed by crazing, i.e. the formation of 

4 interpenetrating micro- voids pp. Similarly, the plastic deformation of semi- 

5 crystalline polymers is strongly correlated to the nucleation of cavities in the 
e amorphous region [2]. Although essential to control deformation and failure 

7 of many organic materials, cavitation in glassy polymers under load is poorly 

8 understood. To our knowledge, the microstructural causes, or the precursors 

9 of cavitation at a microscopic scale, are not clearly identified. Although it is 

10 known that impurities or surface defects aid the nucleation of cavities [21 H] , 

11 it is presently not possible to predict where cavitation will take place in the 

12 polymer. 

13 Classical nucleation theory, where elastic energy is balanced by the cre- 

14 ation of free surface, was used by Argon to model the cavitation nucleation [1] . 

15 Estevez et al. investigated the fracture toughness in glassy polymers using 

16 mechanical approaches with empirical constitutive equations to describe the 

17 competition between shear yielding and crazing [5]. They noted that the 
is development of crazes is favored by a fast local deformation. 

19 According to Gent [6], crazing in glassy plastics can be attributed to a 

20 local stress-activated devitrification. It is generally agreed that large triax- 

21 ial tensile stresses are needed to induce cavitation, which forms the basis 

22 of several macroscopic craze initiation and cavitation criteria [H] • Molec- 
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23 ular dynamics simulations of polymer glasses also found a transition from 

24 shear yielding, which obeys a pressure-modified von Mises yield criterion [9] , 

25 to cavitation as the hydrostatic pressure becomes negative, but have not 

26 yet investigated the connection between cavitation and local microstructural 

27 configuration. More recent simulations explored correlations between the lo- 

28 cation of failure, higher mobility regions and a higher chain ends density [10] 

29 acting then as local defects, or a local, stress-induced disentanglement of 

30 chains [IT]. In the latter work a primitive path algorithm was used to mon- 

31 itor the entanglement network in a sample undergoing triaxial deformation, 

32 and it was found that regions undergoing crazing were also depleted in terms 

33 of entanglements. 

34 The local mechanical properties are a determining factor to understand 

35 the response of systems under strain. Yoshimoto et al. [T2l [13] have cal- 

36 culated the local elastic modulus in a coarse grained polymer glass using 

37 a thermodynamic approach based on stress fluctuation . They found that 

38 polymers are mechanically heterogenous at local scale. Papakostantopoulos 

39 et al. [H] have studied the earliest local plastic events observed in the elastic 

40 regime of polymer glass. They found that these irreversible events take place 

41 in domains that exhibit a low positive elastic modulus. Analogous results 

42 were obtained by Tsamados et al. on Lennard- Jones glasses submitted to a 

43 quasistatic shear strain [TS]. They found a correlation between high nonaffine 

44 displacements and local low elastic modulus. 

45 It is unclear, however, if these criteria can be used in a predictive manner, 

46 in the sense that the cavitation event could be predicted from the configura- 

47 tion of an unstrained system. 
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Figure 1: (a) Stress-strain curve of a glassy polymer at T — 0.2 during a triaxial tensile 
test. A peak in the mean square non-affine displacement is observed simultaneously with 
the drop of the stress due to cavitation. Markers 1, 2, 3 and 4 indicate the selected 
configurations for which distributions of NAD are plotted in (b). These distributions 
show an exponential tail for configuration 4, due to the large amplitude motions caused 
by cavitation, (c) Representation of the NAD in a cross-section of the sample containing 
the cavity, for configuration 4. 

48 In this paper, we will therefore investigate the correlation between the 

49 microstructure of a homopolymer at the segmental level and the nucleation 

50 of cavities, in an attempt to find a microstructural predictor of such events. 

51 Section [2] will present the methodology and will demonstrate that the non- 
52 affine particle displacement (NAD) is a particularly suitable tool for charac- 

53 terizing and locating cavities. Section [3] will be devoted to the investigation 

54 of possible causes of cavitation, namely (i) local excess of free volume, (ii) 

55 local excess of atomic stress, (iii) local density of beads and chain ends, and 

56 (iv) local bulk modulus. We will show that while (i)-(iii) bear little correla- 

57 tion to the NAD, the local bulk modulus (iv) has a much better potential to 

58 predict the cavitation event. 
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59 2. Methods 



60 2.1. Molecular dynamics simulations 

ei Molecular dynamics (MD) simulations were carried out for a well es- 

62 tablished coarse-grained model, in which the polymer is treated as a linear 

63 chain of N beads of mass m, which we refer to as monomers, connected 

64 by stiff anharmonic springs that prevent chain crossing and breaking [16J. 

65 The beads interact through a conventional 6-12 Lennard- Jones potential 
ee that is truncated at 2.5 times the particle diameter (The tiny discontinu- 
67 ity of the force at the cutoff distance, less than 1% of the maximal attractive 
es part, has no consequence). All quantities will be reported in units of the 

69 Lennard- Jones length scale a and energy scale e, and the characteristic time 

70 is Tu = a/hio" 2 Je. Newton's equations of motion are integrated with the 

71 velocity Verlet method and a time step At = 0.006. Periodic simulation cells 

72 of initial size L x (0) = L y (0) = L z (0) = 34.2 containing M = 215 chains of 

73 size N = 200 beads were used with a Nose-Hoover thermostat, i.e. in the 

74 NVT ensemble. All samples were generated using the "radical-like" polymer- 

75 ization method [17J. The polymerization starts from a Lennard- Jones liquid, 

76 where 215 beads are chosen randomly to behave as "radical" sites. Each rad- 

77 ical bead is allowed to connect to a free and nearest neighbor with a strong 

78 covalent bond. The radical sites are then transferred to the new connected 

79 beads, allowing thus the growth of all chains. If no monomers are near the 
so radical, no FENE bond is created. Another attempt will be performed at the 
si next growth stage. Between two growth stages, the entire system is relaxed 

82 during 100 MD steps. The polymerization propagates until all chains reach 

83 their target length of 200 beads. When the generation is terminated, residual 
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84 single beads are removed and the system is relaxed for 10 7 MD steps in NPT 

85 ensemble at T = 1 and P = to reach an equilibrium state. The equilibra- 
se tion leads to a "mean square internal distance" very close to the function 
87 given by Auhl et al. [IB]. The polymer is then rapidly quenched into the 
as glassy state at a temperature T = 0.2 in NPT ensemble (cooling rate: le 

89 per 10 6 MD steps). The glass transition temperature is T g = 0.43 Q The 

90 pressure remains zero and the sample density reaches 1.04 before applying 

91 the deformation. 

92 Triaxial tensile test conditions were employed [19J. The samples were 

93 subjected to a sequence of deformation-relaxation steps, composed of (i) a 

94 rescaling of the simulation box in the tensile direction (y in our case so that 

95 the true strain e yy (t) = ln(L y (t) / L y (0)) , whereas the two other dimensions 

96 remain unchanged, followed by (ii) an MD step in the NVT ensemble. The 

97 deformation rate was chosen to be L y = 0.0025, so that the initial strain 

98 rate is e yy (0) = 7.3 x 10 -5 . Over the range of strain investigated, the true 

99 strain rate remains essentially constant. Note that the applied deformation 

100 trajectory leads to a high level of triaxiality, which is the basic ingredient for 

101 cavitation. As the deformation proceeds, configurations were recorded along 

102 the trajectory in order to analyze their microstructure. 

1 T g has been determined by the slope change observed when the sample volume is 
plotted with respect to the temperature during cooling from T = 1 to T = 0.0001 under 
the NPT ensemble 
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103 2.2. Non-affine displacement: a tool for characterizing cavitation 

104 The mechanical behavior of our glassy polymers under triaxial tensile 

105 conditions is illustrated in Figure [T|a). Three main regimes can be distin- 
loe guished: (i) elastic, (ii) viscoelastic, and (ii) drawing regime, which occurs 
107 at constant stress. In the elastic regime, the increase of deformation will 
loa slightly shift the bead positions from their local energy minima, resulting 
109 in reversible behavior. This regime is limited to a very low strain 0.001 as 
no demonstrated by Schnell [20]. In the viscoelastic regime, stress is relaxed by 
in inter-chain sliding. This stage is limited by a strong drop of stress. When 
n2 a critical deformation is reached, cavities will nucleate and then part of the 
in stored elastic energy is released as free surfaces open up. Note that the strain 
H4 hardening regime is not shown in Figure [TJa), since it occurs at larger strains 
H5 when the entanglement network of chains and fibrils becomes stretched [9]. 
lie The detection of cavity nucleation could be performed visually on snap- 
H7 shots that are regularly stored during the course of the tensile test. How- 
n8 ever, small cavities in a three dimensional sample can be delicate to observe. 
n9 Therefore, a more versatile indicator is needed. The non-affine displacement 

120 (NAD) is the perfect candidate for such observation and has been success- 

121 fully used to monitor local plastic activity in 2D amorphous Lennard- Jones 

122 packings under athermal quasistatic deformation [21J. Note that NAD fluc- 

123 tuations can not find their origin in the thermal motion of atoms since, in 

124 the framework of this paper, specimens are maintained well below their glass 

125 transition temperature (T = 0.2 < T g = 0.43). Moreover, the NAD can be 

126 used as a routine tool and it starts to increase locally, in the early stages 

127 of cavity nucleation, even before the cavity could be observed visually on a 
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128 snapshot of the sample. 

129 The non-affine displacement (u^ a ) is defined as the difference between 

130 the mean displacement of a bead % during time St (r*(t + St) — and 

131 the mean displacement it would experience if the deformation were perfectly 

132 affine, i.e homogeneous at all scales, 

uUt) = r l (t + St) - r\t) - e yy {t)St rj(f)e„, (1) 

133 where r l {t) is the position of bead i at time t, r y (t) is the projection of 

134 this position along along the y axis and St is the time elapsed between two 

135 configurations where the NAD is evaluated (typically 30r). 

136 In Figure pja), the cavity nucleates at e yy = 0.09. At the same strain, 

137 the NAD exhibits a peak. Beads that exhibit the largest NAD are those 

138 which belong to the surface of the cavity (see Figure |l[c)). Figure [IJb) 

139 shows the evolution of the NAD distribution for several deformations. Before 

140 cavitation, increasing the deformation shifts the distribution tail to larger 

141 NAD until the very moment at which the growth of a cavity occurs, which is 

142 associated with very large values of NAD (see Figure [T](c)). A threshold for 

143 NAD magnitude has been defined: if |u^J > 1.5a at the yielding point, the 

144 bead i is said to belong to the cavity surface. This threshold is used to identify 

145 the "cavity beads" in order to follow some of their local properties. The 
we position of the cavity is defined as the centre of mass of these "cavity beads" . 

147 After cavitation, the distribution returns to a narrower shape. Note that the 

148 NAD distribution broadens even before the stress drop in the stress-strain 

149 curve, due to the nucleation of the cavity. In the following sections, NAD 

150 will be used as a quantitative tool for investigating the possible correlations 

151 with other microstructural or mechanical properties, such as Voronoi volume, 
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Figure 2: Voronoi volume distributions of configurations extracted over the deformation 
trajectory. Cavitation can be clearly seen in the tail of distribution. After cavitation, the 
distribution reverts to a narrower shape. 

152 hydrostatic stress, local density and local moduli. 



153 3. Microstructural causes and precursors of cavitation 

154 In this section, we will attempt to correlate NAD fluctuations with some 

155 local properties measured at the scale of a single "atom" ( Voronoi volume 

156 and stress per atom), and properties averaged on the scale of a few particle 

157 diameters (chain end density and bulk modulus). 
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158 3.1. Voronoi volume fluctuations 

159 The concept of free volume has been extensively used to explain many 
wo specific properties of supercooled liquids and glasses. Free volume is defined 
lei as the volume in excess compared to an ideal disordered atomic configuration 

162 of maximum density. One of the simplest way to compute free volume on a 

163 local scale (and to avoid the ambiguity of the above definition) is the Voronoi 

164 tessellation, which uniquely assigns a polygonal volume to each bead, formed 

165 by intersecting the planes bisecting the lines between different bead centres, 
lee In order to determine whether local fluctuations of free volume (or Voronoi 
167 volume) favour the nucleation of a cavity, we used the voro++ routine to 
lea calculate the volume associated to each bead 

169 Voronoi volume and deformation level. . Figure [2] shows the effect of the de- 

170 formation on the Voronoi volume distribution. Increasing the deformation 

171 will increase almost homogeneously the free volume until cavitation takes 

172 place. During and after cavitation, the Voronoi volume distribution exhibits 

173 a significant tail representing the beads belonging to cavity walls. Note that 

174 after cavitation, the distribution relaxes to a narrower shape. Therefore, the 

175 cavitation process can be seen as an event, which localizes or precipitates the 

176 excess of free volume introduced by deformation. 

177 Voronoi volume and beads functionality. . Figure [3] compares the mean Voronoi 

178 volume evolution of both regular beads and chain ends. It can be seen that 

179 chain ends exhibit a larger Voronoi volume, which is not surprising since, by 



See http : //math . lbl . gov/voro++/| and ref . [22] > where a very early version of this 



code was used. 
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Figure 3: Evolution of the mean Voronoi volume of chain ends and other beads. Inset: 
Voronoi volume distributions during the course of deformation. No correlation can be 
observed between chain ends and Voronoi volume variation during tensile test. 
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180 construction, covalent and Lennard- Jones bonds have their energy minimum 

181 at 0.9a and 1.12a, respectively. Note that when cavitation occurs, the mean 

182 Voronoi volume of chain ends becomes very noisy due to statistical limita- 

183 tions. The insets of Figure [3] show that the Voronoi volume distributions 

184 have a Gaussian shape, which shows the presence of low Voronoi volumes 

185 (much lower than the volume of an ideal disordered configuration). This 
we calls into question the very concept of free volume, which is defined as that 
is? part of the atomic volume that can be redistributed throughout the system 
las without change in energy [2U [25], i.e. the volume of an ideal disordered 
189 configuration. These points of extremely low volume could be related to the 
wo constriction points introduced by Stachurski [23J and, to a larger extent, to 
hi the quasi-point defects of Perez [TJ , which represent points of high fluctuation 

192 of free energy. 

193 Voronoi volume and non-affine displacement.. Within the free volume ap- 

194 proach, deformation induced relaxations are supposed to be correlated with 

195 the available free volume. Zones of larger free volume will therefore deform, 
we changing the potential energy landscape and providing more free volume to 

197 zones of initially larger free volume. This explanation is often proposed to 

198 describe the formation of mechanical instabilities such as cavitation or shear 

199 bands. Motivated by these ideas, the search for a relationship between the 

200 magnitude of the NAD and the Voronoi volume becomes relevant. 

201 Figure [4] shows a scatter plot obtained during deformation, where the 

202 magnitude of the NAD and the Voronoi volume were taken as variables. This 

203 scatter plot does not show a clear tendency for a correlation between NAD 

204 and Voronoi volume. Free volume represents the potential space for motion 
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Figure 4: (Color online) NAD magnitude of beads plotted against their Voronoi volume at 
several strains e=0, 0.03, 0.06 and 0.09 (cavitation). Arrows describe the higher density 
value of each variable. No correlation was found between these variables for all strains. 
"Cavity cluster beads" are marked by the symbol o: no remarkable trend can be dis- 
tinguished, except during cavitation, where they exhibit larger NAD and slightly larger 
Voronoi volume. 
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205 but it can not be seen as being a causal factor of the NAD and cavitation. The 

206 "cavity cluster" beads are also shown in this plot. In both cases, the points 

207 are distributed randomly and no noticeable trend was found, except during 

208 cavitation, where these beads exhibit larger NAD and slightly larger Voronoi 

209 volume. This analysis (not shown in this paper) was performed for several 

210 other temperatures (T = 0.01 and T = 0.1) and no correlation was found 

211 under these conditions either. Note that before cavitation, the magnitude of 

212 NAD remains much less than inter-atomic distance (a), in other words, the 

213 deformation is purely affine. 

214 3. 2. Stress fluctuations 

215 The local stress on any given bead can be obtained by dividing the classi- 

216 cal expression of the virial stress by the Voronoi volume V m of atom m [26J , 



T • r ''"' -17) ■ ( 2 ) 



218 where v™ is the velocity ith component of atom m; / J mn and r™ n are 

219 the zth component of force and distance between two interacting atoms m 

220 and n, respectively. The first term of this equation represent the kinetic 

221 contribution and the second one is the Cauchy stress. The hydrostatic stress 

222 Shyd was calculated by computing the trace of the stress tensor, Shyd = 

223 -((Tn + cr 22 + 0"33)/3 

224 Figure [5] compares the distributions of hydrostatic stresses at several 

225 strains during deformation. In the initial undeformed configuration, the dis- 

226 tribution shows an exponential tail towards negative values. As the deforma- 

227 tion increases, the negative values of the hydrostatic stress are progressively 
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Figure 5: Distribution of hydrostatic stress at different deformation levels. The distribu- 
tion narrows until the cavity opens, then broadens again. 

228 relaxed, so that the distribution narrows and becomes more symmetrical 

229 just before cavitation takes place. These results are consistent with those 

230 of ref [2S]- After the cavitation has occurred, large negative values of the 

231 stress are again obtained, since the average free space is decreased as shown 

232 previously in Figure [2} 

233 In order to investigate the correlation between NAD and the atomic hy- 

234 drostatic stress, a scatter plot of these quantities is displayed in Figure [6] 

235 Apparently, there is no direct trend for a correlation between NAD and the 

236 hydrostatic stress at the scale of individual beads. The hydrostatic stress was 
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237 also evaluated by considering each contribution separately (pair, bonded) at 

238 several strains, and again no correlations were found. When specific beads 

239 (chain ends and cavity cluster) values are selected in these scatter plots, the 

240 corresponding points appear to be a randomly chosen subset of the total 

241 sample. This absence of correlation may appear surprising, as the presence 

242 of a high local stress is often expected to result in plastic deformation. Note 

243 however that this result is consistent with a recent study [2S] which showed 

244 no correlation between atomic stresses and shear yielding in polymers. In 

245 an analogous way, a previous study on sheared glasses [27] also failed to 

246 find a direct correlation between local stresses and the relevant local plastic 

247 deformation (shear transformations in that case). 

248 It may be, however, that a more coarse grained characterization is nec- 

249 essary to identify such correlations, and that the cavitation events are the 

250 result of a local heterogeneity that extends beyond the scale of individual 

251 beads. In order to assess this hypothesis, we describe briefly in the next 

252 section studies performed on density fields defined at a more coarse grained 

253 scale. 

254 3.3. Coarse grained densities 

255 The opening of a cavity under strain can be seen as a collective event, that 

256 involves at least those atoms that will form the cavity "skin" at the end of 

257 the process. The corresponding mechanical instability may therefore be the 

258 result of some density anomaly that extends over a region larger than a single 

259 atom size or Voronoi cell. We therefore have also explored the properties of 

260 our polymer system on such a coarse grained scale by defining continuous 

261 fields from the atomic positions. Various possibilities are available for such a 
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Figure 6: (Color online) Scatter plot of the non affine displacement against hydrostatic 
stress shown at the same deformation as in Figure [4] No clear trend for a correlation can be 
established. The values corresponding to the beads that surround the cavity labeled by the 
symbol o are randomly dispersed, thus preventing one to identify any specific correlation 
for these beads. 
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Figure 7: (Color online) Evolution of mean local density calculated within cubic bins of 
size 5 x 5 x 5 cr 3 . The average of the overall voxels (dashed line) is compared with the 
average of voxels that are associated with the cavity cluster beads (solid line). Symbols (-) 
denote the upper and lower values for each curve. The spatial distribution of local density 
taken at e = 0.08 is shown in (b) and the corresponding map of chain end density in (c). 
The arrows describe the non affine displacement of the "cavity cluster beads" . The dark 
spots identify the low density vicinities (b) and the high chain end density (c). 

262 coarse graining procedure [2Hl US! EH] ; here we choose the simplest one, which 

263 consists in computing the densities on a regular grid by assigning to each grid 

264 node the atoms that belong to a fixed "voxel" volume around this node. The 

265 voxel size is taken in the range 5a to 7a, which was shown in similar studies 

266 [301 ESI EH] to permit a good description in terms of continuous fields (with 

267 about 120 monomers per voxel) while preserving the locality and possible 

268 spatial heterogeneity of the variables under consideration. 

269 We have attempted to coarse grain and to correlate with the appearance of 

270 cavities two of the densities examined previously at the atomic level, namely 

271 the density of chain ends and the density of monomers. The local density 

272 field is defined as: 

Pi = y^V ( 3 ) 

273 where rij is the number of beads within a voxel i, and is the Voronoi 
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274 volume of the bead j that is included in the voxel i. Figure u[a) shows the 

275 evolution of this local density with strain, and compares the average value 

276 with the value observed in the vicinity of the cavity. 

277 These data show that the local density in the vicinity of the cavity fol- 

278 lows the mean value until a deformation of e = 0.06. Although cavitation 

279 does not occur until e = 0.09, the local volume begins to decrease earlier 

280 (see inset). This observed trend can be interpreted by the fact that cavi- 

281 tation starts earlier than the drop of stress in the stress-strain-curve. This 

282 "pre-cavitation" behavior can be interpreted as resulting from a dynamical 

283 equilibrium between the elastic energy and the free surface energy of cavity 

284 with relatively small radius. This situation remains stable until the cavity 

285 reaches a critical radius (roughly estimated 2a), beyond which the size of 

286 the cavity increases rapidly. The spatial distribution of the local density at 

287 e = 0.08, just before the opening the cavity is shown in part (b) of Figure 

288 [7j The lowest density spots are far from the expected position of the cavity, 

289 but a low density can be noticed in the cavity vicinity. In general, we have 

290 checked that a systematic decrease in density prior to cavitation is specific of 

291 the points that are located in the vicinity of the emerging cavity. Other points 

292 may display fluctuations in their values of the density, but these fluctuations 

293 remain uncorrelated with cavitation events. After the cavity nucleation, the 

294 low density regions that did not form cavities release their excess free volume 

295 introduced by the triaxial deformation condition. Therefore the local density 

296 of such regions return to values similar to regions that are not involved in 

297 the cavitation. In conclusion, local loss of density should be seen rather as a 

298 consequence than as a cause of cavity nucleation. 
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299 As was mentioned above (in section 3.1 ), the free volume was found to be 

300 correlated with the bead connectivity. Chain ends exhibit a higher Voronoi 

301 volume compared to other monomers, and a lower density of beads could be 

302 expected where a higher density of chain ends is present. We therefore define 

303 a local density of chain ends p°- E - as 

y~.CE. 

P?- K = \, (4) 

304 where nf' E ' is the number of chain ends within a voxel i, and V{ is the volume 

305 of the voxel. Figure [7](c) shows that, at this level of coarse graining, the spa- 

306 tial distribution of chain ends is uncorrelated with the local density of beads 

307 and also with the cavity position. This indicates that the modification of the 

308 packing density by the presence chain ends is insignificant. Summarizing, the 

309 coarse grained density of beads exhibits a limited success as a predictor for 

310 cavity formation, as its evolution can be correlated with the formation of a 
3n cavity only shortly before the event actually takes place. The coarse grained 

312 density of chain ends, on the other hand, does not correlate well with the 

313 total density or with cavitation. 

314 3.4- Local mechanical properties 

315 Our last attempt to identify a microstructural predictor for cavitation 

316 events is inspired by previous work on simple glassy systems under shear 

317 deformation, in which a low value of the shear moduli was identified as a 

318 good indicator for the occurrence of the relevant local plastic events, shear 

319 transformation zones [151 EH] • Here the relevant events involve a local dilata- 

320 tion of the material which eventually gives rise to a cavity, and points to the 

321 local bulk modulus possible predictor. 
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322 Local heterogeneity in the elastic properties of glasses is now a well doc- 

323 umented feature, with a number of studies having shown that the moduli 

324 defined at intermediate scales (of the order of 10 atomic sizes) are those of 

325 an isotropic but heterogeneous material. At such scales, a typical glassy 

326 sample can be described as a consisting of coexisting "hard" and "soft" re- 

327 gions. This behavior is independent of the precise method which is used to 

328 define the coarse grained elastic constant, which may involve either the use 

329 of statistical mechanical formulae at a local scale [131 E] , or exploiting the 

330 linear relationship between coarse grained stress and strain field [15] . Here 

331 we present results for the local bulk modulus obtained from a third approach, 

332 originally introduced by P. Sollich et. al. [ST] , which has the advantage of be- 

333 ing easily implemented at a reduced computational cost. The method can be 

334 summarized as follows: one first defines a coarse graining volume as a Active 

335 shape that encapsulates a number of beads. The shape was chosen spherical 

336 in order to reduce any potential boundary effects, and the radius equal to 

337 3.5 particle diameters, consistent with the typical coarse graining scales used 

338 in other methods [T5l [T3J. The entire sample is then deformed affinely (in 

339 this case using a uniform dilation of all bead coordinates). After this ho- 

340 mogeneous deformation, all beads are kept frozen, except those contained in 

341 the coarse graining volume which are allowed to fluctuate in a constant vol- 

342 ume, constant temperature molecular dynamics trajectory (here we perform 

343 a trajectory at a rather low temperature, T = 0.01). The increase in the 

344 hydrostatic stress S™ yd within the spherical volume m is obtained from the 

345 virial stress formula, and the local modulus K m can be defined by dividing 
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Figure 8: Density distribution of local bulk modulus for the same specimen at several 
deformation levels. Distributions arc shifted by their mean values. The inset shows the 
evolution of the mean values. Error bars are deduced from the standard deviation. 

346 this stress by the imposed increase in volumetric strain $ m : 

*- = 'W- <«> 

347 In order to improve the accuracy on K m , it has then been averaged over 

348 a dozen expansion tests within the domain < d m < 10~ 5 . A sequence of 

349 deformation (isotropic expansion) and relaxation steps is applied over the 

350 sample, the gauge volume V m and S™ yd are measured after each step. The 

351 expansion is limited to a very low deformation amount since the measurement 

352 is restricted in the elastic regime only. Substituting the i? m by its definition 
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353 



dV m 
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leads to another form of equation (p 
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dV 



(6) 
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354 or equivalently 



K m -\n{V m )\l = S] 



im 1 1 
frj/dlO' 



(7) 



355 where V m and S™ yd are integrated along the deformation trajectory from 

356 to t. This method allows us to obtain an accurate determination of K m 

357 in the linear regime by fitting the data obtained for S , ^ d =f(Zn(l /m )). This 

358 procedure was applied along each tensile deformation trajectory, for positions 

359 of the center of the coarse graining volume distributed on a regular grid. 

360 Figure [8] compares the statistical distributions of the local bulk modu- 

361 lus at different strain levels. The plotted distributions are shifted by their 

362 mean values to facilitate comparison of their shapes. Curves remain sym- 

363 metrical and Gaussian, whatever the applied strain before cavitation. As the 

364 deformation increases, the distribution will become slightly narrower. This 

365 behaviour is consistent with the trend described in the previous sections, that 

366 the polymeric system tends to homogenize its local stress under an applied 

367 deformation. The mean value of the local bulk modulus (see inset) decreases 

368 continuously as the deformation increases and more free volume is introduced 

369 in the system. 

370 In order to investigate the spatial distribution of the local bulk modulus, 

371 two-dimensional slices in a plane perpendicular to the tensile direction were 

372 taken at the level at which the cavity is observed. Figure [9] shows a sequence 

373 of such bulk modulus cartographies that are captured along the deformation 

374 trajectory. Each slice corresponds to one of the blue markers on the stress 
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Figure 9: (Color online) XY Layers taken from the spatial distribution of the bulk modulus 
(tensile direction is Z in this case). The Z coordinate of all layers corresponds to the 
cavity position. Every layer is indexed (bottom left corner) referring to its position in the 
stress-strain curve on the upper left. Theses maps are superimposed with the nonafHne 
displacement ate = 0.09 that shows the formation of the cavity. 
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375 strain curve (first plot in figure [9]). The nonaffine vectors describing the 

376 formation of a cavity are also plotted on each map. As can be seen, the local 

377 bulk modulus fluctuates between high and low values at each strain, and 

378 the position at which the cavity appears corresponds to one of the low bulk 

379 modulus sites identified in the starting configuration. When the deformation 

380 increases, an extremely low value of bulk modulus appears in the expected 

381 position of cavity, as in slices (4) and (5). The lower value indicated here is 

382 not only a local minimum in the plane of the figure, but instead corresponds 

383 to the lowest value for the entire sample. 

384 In the light of this strong correlation between NAD and elastic modulus, 

385 the cavitation process in glassy polymers can be described in the following 

386 manner: The polymeric system exhibits some fluctuations in the local elastic 

387 bulk modulus. As the deformation progresses, the statistical distribution of 

388 the bulk modulus changes: The mean value decreases, but the contrast of 

389 spatial distribution is conserved. At relatively high strain, one of the zones 

390 that initially displayed a low bulk modulus will reach an anomalous value, 

391 resulting in a favorable location for the subsequent growth of a cavity. 

392 We will now investigate whether this behavior should be described as 

393 deterministic (the cavity systematically forms in a particular zone) or rather 

394 statistical (the cavity forms randomly in one of the zones with a low modulus) 

395 process. To this end, the same system was subjected to three tensile tests 

396 with different tensile directions x, y and z. The positions at which cavitation 

397 takes place were recorded and compared after each test. 



399 cleate in different zones. The same behaviour was also found for several 



398 




shows that, for the same initial configuration, cavities nu- 
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400 systems with different temperatures. The cavities systematically nucleate in 

401 zones that are characterized by a low modulus in the initial state, however 

402 the specific site at which it is observed depends on the deformation path and 

403 on the tensile direction. 

404 4. Conclusions 

405 In this work the relationship between a cavitation event in a glassy poly- 

406 mer undergoing a tensile test and the local properties has been investigated 

407 with molecular dynamics simulations. Several properties have been analyzed 

408 at two different length scales: the elementary scale of the monomer, and a 

409 coarse graining scale of 5 to 10 particle diameters. Independent of the scale 

410 under consideration, we find that the density of monomers or the density of 

411 chain ends do not correlate with the subsequent appearance of a cavity. In 

412 contrast, the bulk modulus in the unstrained configuration displays fluctu- 

413 ations that can be directly related, in a statistical sense, to the appearance 

414 of a cavity at large deformations. Note that very similar conclusions were 

415 reached by Toepperwein and de Pablo in a recent study that considered both 

416 homopolymers and composite systems [32]. This situation resembles those 

417 observed in glassy materials under volume conserving shear, where a weak 

418 shear modulus indicates a tendency for plastic rearrangement. 

419 
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Figure 10: (Color online), (a) Positions of cavities in a single sample after tensile tests 
with different straining directions. The vectors describing the formation of a cavity are 
superimposed with the corresponding map of local elastic bulk modulus imaged in the 
initial, undcformed, configuration. Slices (b), (c) and (d) are located in the corresponding 
position of the cavity nucleation and they are perpendicular to the tensile directions X, Y 
and Z, respectively. 
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